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ABSTRACT 

The conventional method of generating initial conditions for cosmological A^-body simulations in¬ 
troduces a significant error in the real-space statistical properties of the particles. More specifically, 
the finite box size leads to a significant underestimate of erg, the correlation function, and nonlinear 
effects. I implement a method of generating initial conditions for A-body simulations that accurately 
models the real-space statistical properties, such as the mass variance in spheres and the correlation 
function. The method requires running ensembles of simulations because the power in the DC mode 
is no longer assumed to be zero. For moderately sized boxes, I demonstrate that the new method 
corrects the underestimate in the mass variance in spheres and the shape of the correlation function. 

I also argue that subtracting Poisson noise from the power spectrum is a dangerous practice. Code to 
generate initial conditions to second order in Lagrangian perturbation theory with the new method is 
available at http://www.astro.princeton.edu/~esirko/ic. 

Subject headings: cosmology: theory — large-scale structure of universe — methods: A-body simula¬ 
tions — methods: numerical 


1. INTRODUCTION 

Cosmological A-body simulations are an important 
tool in understand ing nonlinear structure formation 
(|Rertsch in srerl . For a consistent framework, this 

paper focuses on the collisionless variety, where each par¬ 
ticle represents a sampling of the underlying dark matter 
collissionless fluid and baryon physics is assumed to be 
negligible on large enough scales. However, the main re¬ 
sults of this paper can be applied to any type of numerical 
simulation, such as hydrodynamic simulations of galaxy 
formation and Lya clouds, in which the mean density 
is representative of a cosmologically relevant part of the 
universe. As with any integrator, the generation of ini¬ 
tial conditions is as important as the integration itself. 
However, there seems to be a relative dearth of informa¬ 
tion on initial conditions to A-body simulations in the 
literature, besides the standard “we use the Zeldovich 
approximation” tagline. This paper is an attempt to fill 
that gap. 

The app roach for the present work was motivated by 
Pen pointed out that in order to accurately 
model certain statistics — those most directly connected 
to analytic models of nonlinear structure formation — 
one should ensure that the real-space statistical proper¬ 
ties (e.g., the correlation function) of the initial condi¬ 
tions are accurately modeled, instead of the fc-space sta¬ 
tistical properties (e.g., the power spectrum) as is usually 
done. In this paper, I denote initial conditions gener¬ 
ated by the former method ^-sampled IC’s (because they 
accurately model the correlation function), and by the 
latter, P-sampled IC’s (because they accurately model 
the power spectrum). In the limit that the box size 
L —f oo, the two methods converge to each other, but 
in general, it is impossible to accurately model the corre¬ 
lation func tion and power spec trum simultaneously in a 
finite box. iBertschiugerl (|2flfllH released a multiscale ini- 
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tial conditions g enera t or cal led GRAFIC2. While building 
upon ideas from iPeiJ l|1997ll . GRAFIC2 does not actually 
incorporate Pen’s crucial suggestion of using ^-sampled 
IC’s. Bertschinger claimed that using ^-sampled IC’s 
would be inconsistent with A-body codes that use peri¬ 
odic boundary conditions, but there is actually no such 
inconsistency because P- and ^-sampled IC’s share the 
properties of periodic boundary conditions and discrete 
wavenumbers. 

Historically, simulations were run until the mass vari¬ 
ance in spheres erg reached a certain value, and that 
epoch was called “today.” But with the advent of “preci¬ 
sion cosmology” from accurate CMB measurements, tig is 
no longer a fudge factor; it is a well-defined cosmological 
parameter that should be input into the initial conditions 
as a starting point to simulations. Speaking broadly, the 
advent of precision cosmology requires simulations to be 
accurate at the percent level. The motivation for this 
paper, therefore, is to get A-body simulations accurate 
enough for precision cosmology. 

The discrepancy between P- and ^-sampled IC’s is 
closely related to the problem of finite box size. By 
studyi ng the fraction of mass in halos, iBagfa Ki. Ravi 
il2?inl have recently argued that box sizes should be 
much larger than those typically used. But their analysis 
is for P-sampled IC’s. While I do not use the statistic of 
halo mass fraction in the present study, the analysis of CTg 
in m below implies that using ^-sampled IC’s solves the 
problem of underestimated real-space statistics in finite 
b oxes “to first order.” 

I.Iovce. Levesque, fc Marco^ <|2004jl have proposed a 
new method of generating initial conditions by modify¬ 
ing the potential of a one-component plasma and inte¬ 
grating trajectories of particles (similar to generation of 
glass initial conditions). They claim that this method 
will accurately model the correlation function and power 
spectrum simultaneously. A clarification of this apparent 
contradiction is in order. Their paper is an attempt to 
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solve an orthogonal problem, that of ignoring finite box 
size effects but focusing on nonzero Ax = LjN^I'^ ef¬ 
fects. In fact, particle discreteness has two adverse conse¬ 
quences. Fi rst, nearest-neighbor effects are dynamicall y 
important llBaertschiger. Joyce, fc Svlos Labinl I2002D . 
This consequence of particle discreteness seems impos¬ 
sible to alleviate without increasing the mass resolu¬ 
tion of a given simulation. Second, the correlation 
function of a perturbed lattice of particles is highly 
oscillatory and drastically misrepresents the true cor- 


relation function of the underlying continuous fluid 

j|Bag-tediij,cr_&^jyjns Ljahmi20(^ fc Domingued 

I2003H, Joyce fc Marcosll2004|l . The present paper ignores 

particle discreteness effects, instead focusing on the prob¬ 
lem of finite box size L. However, as seen in below, 
the oscillations in the correlation function for a lattice are 
erased cleanly as the simulation evolves and the under¬ 
lying lattice disappears. It is reasonable to suppose that 
the effects of the oscillatory corre lation function (which 


l.lovce. Levesone. fc Marcoa l|2004ll attempt to solve) are 
not as important as the effects of a systematically under¬ 
estimated correlation function (which the present paper 
attempts to solve). 

As we will see, .^-sampled IC’s introduce a nonzero 
variance of the DC mode of a simulation. This im¬ 
plies that running an ensemble of simulations will be 
necessary, where each realization’s DC mode is drawn 
from an appropriate distribution. Running an en¬ 
semble is a good idea anyway, because of the small 
number of independent modes (~ 3) with wavelength 
~ L. To combat this l a tter e f fect of sample vari¬ 
ance. iKnebe fc D omingu^ ll2003ll. iBaugh fc Efsta thioiJ 
lll99'^ . and lBaiigh!*Caztan^^fc EfstathiouTilTQQiTfl ran 
ensembles of 10 simulations but used P-sampled IC’s, 
ignoring the DC mode. (The latter authors compen- 
sated for the lack of DC mod e by using the technique of 
iCouchman fc CarlbergI lll992tl of extrapolating the power 
spectrum at small k by linear theory.) In contrast, var¬ 
ious authors have proposed methods of extracting an 
“ensemble’s worth” of information out of a single re¬ 
alization. For example, a realization, being periodic, 
can be replicated and longer-wavelength inodes linearly 
addedjTormm fc Bertscliin geilll99fiHColjll997ll . Simi¬ 
larly, [Couchma^^^CaHbe^'^^filpInked nonlinear evo¬ 
lution at high k with linear theory predictions at small 
k in their simulations. While these methods are useful 
in their own right, in this paper I will assume that com¬ 
putational r esources are availab le to run an ensemble of 
simulations. iFrenk et al.l ill 988fl ran an ensemble of sim¬ 
ulations very similar in spirit to the methods outlined in 
this paper, using P-sampled IC’s for k ^ 0, but taking 
into account the DC mode. The present paper extends 
their method in several ways: the DC mode of each re¬ 
alization should be drawn from a Gaussian distribution, 
and all modes of the power spectrum, not just the DC 
mode, should be convolved with the window function for 
^-sampled IC’s. _ 

The next subsection, o is an exposition to some 
of the terminology and definitions used in this paper. 
I outline the conventional method for generating initial 
conditions with the Zeldovich approximation in m SI 
discusses the problem with the conventional method, and 
(SJdiscusses the solution of using ^-sampled IC’s. fJHcon- 


siders the DC mode, an important complication arising 
from ^-sampled IC ’s. m presents numerical simulations 
designed to test .^-sampled IC’s. §(0151 analyze certain 
statistics of the particle distributions: the mass variance 
in 8 /i“^Mpc spheres, the correlation function, and power 
spectrum, respectively. A comment on subtracting the 
Poisson noise contribution follows in HIOI illll is the sum¬ 
mary. 


1.1. Definitions, notation, preliminaries, etc. 

As is usual, we will assume that the cosmological sim¬ 
ulation tracks the evolution of a periodic comoving cu¬ 
bical box with side length L, so that the simulated uni¬ 
verse is at least homogenous out to infinity. The vol¬ 
ume V = L^. There are N collisionless particles in the 
box, representing a sampling of the dark matte r field . 
The early cosmological simulations of iDavis et al IdUll) 
used N = 32^, and present-day state-of-the-art simula- 
tions by the Virgo co nsortium have used 2160^ particles 
llSnringel et al.ll200fiD . 

The density field p, whether continuous or discrete, 
becomes S = pjpo — 1 in dimensionless form, where po 
is the average density. If the density field consists of 
point particles, (5(x) = ^ ^(5 d(x — Xp) — 1 where (5 d 
is the Dirac delta function and Xp is the location of 
each particle. The Fourier t ransform equations are (see 
iHocknev fc EastwoodI lll988li . Appendix A) 

J(k) = J <5(x)e-**^'"d3x (1-a) 

<5(x) = i^5(k)e**^- (1-b) 

Note equation OI can be derived by approximating the 
Fourier transform integral (5(x) = j2^ I 'j(k)e®’^ ’*’d^k. 
A^-body codes follow the trajectories of particles in a La- 
grangian framework. Despite the discretization of matter 
into particles, the density field is continuous. But it oc¬ 
cupies a finite volume, so the integral in e quati on ll-al is 
over the box volume, the sum in equation ll-bl is an in¬ 
finite series, and k is forced to have values S {i,j,k)Ak 
where Afc = 27r/L and (i,j,k) are integers G (— 00 , 00 ). 
For comparison, a Eulerian discretization would turn the 
inte gral i n eauation ll-al into a sum, and the sum in equa¬ 
tion ll-bl would be Nyquist limited. The interpretation 
of the density field given by equation ^ is “literal,” in 
that the particles are delta functions, and it is the con¬ 
vention I will use throughout this paper, but it is worth 
keeping in mind two related issues. First, A-body codes 
use force softening to ensure that the evolution is col¬ 
lisionless. Force softening imparts an effective particle 
shape, so this can be modeled as convolving the “lit¬ 
eral” density field wit h the particle shape. Thus the 
Fourier series leauation ll-bll will effectively be limited at 
large wavenumber, but not Nyquist limited. Force soft¬ 
ening, being a modification to the inverse-square force 
law, is often modeled as an e ffective particle shape in¬ 
teracting with a point mass l|Atha,na,ssonla et al .1120091: 
iHernauist fc Kat3ll989|l . However, in order to explore 
the effects of force softening on the power spectrum, an 
effective particle shape must be found that would be 
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consistent with the modified force law if two of these 
nonzero-sized particles were interacting. The issue of 
force softening’s effect on the power spectrum will not 
be considered further in this paper on the grounds that 
it is unlikely to shed much light on the regions of interest 
of the power spectrum. Second, discretized A^-body fields 
(for cosmological collisionless simulations) are always in¬ 
terpreted as sampling the underlying density field, which 
makes information about the power spectrum at wave¬ 
lengths less than the typical particle separation meaning¬ 
less. However, it remains unclear if there is any “best” 
way to make the conversion from a particle field to a 
continuous one, besides smoothing with a large enough 
filter. In any case, highly nonlinear effects at very high k 
are expected not to have much effect on t he evolution of 
the modes of interest in a simulation (e.g., iPeeble^llflSnl 
§28). 

The power spectrum of the finite-size box is H(k) = 
y(|5(k)p) where the brackets indicate an ensemble av¬ 
erage. Since the universe is isotropic, P{k) = P(k). I 
define the isotropic power spectrum of a specific realiza¬ 
tion of a set of particles in a box with the same formula 
except with the average being over |k| = k, although 
there will be sample variance because of the finite num¬ 
ber of modes. Note that the power spectrum has units of 
volume, because it is a power spectral density in A:-space. 
iBertschinged lll992ll has emphasized that this convention 
should be used. 


2. THE CONVENTIONAL METHOD FOR GENERATING 
INITIAL CONDITIONS 

The goal of an initial conditions generator for a cos¬ 
mological A^-body simulation is to faithfully reproduce 
the statistical properties of the density field in the early 
universe with a finite number of point particles. The 
interpretation of a finite set of point particles as an 
underlying continuous density field is somewhat ambigu¬ 
ous; hence, there are many ways to make this conversion 
and it is not clear if th e re is a ny “correct” answer 
(e.g., iHernquist fc Katzl 119891: iBertschinger fc 
1^ iWeinberg. Hernauist. fc f^; 

Pehigess^j_Schaa£j&vande_Wexgaer^ 1200,1 


Ascasibar fc Binnevll200,^ . We have thus come across 


the first thorny issue of initial conditions generation, 
before even discussing cosmology: how does one convert 
a uniform density field into a set of discrete particles? 
The resultant particle dis tribution can appropri ately be 
called “pre-initial” (e.g., I.Tovce fc MarcosII2004H . It is 
simplest to use a lattice of particles, where the position 
of each particle is given by [{i,j,k) 0.5]L/A^^/^; 

{i,j,k) are integer indices G [0,Af^/^). One could also 
use glass pre-initial conditions, obtained by evolving 
randomly placed particles with the sign o f gravity re¬ 
versed (with damping to help convergen ce) llWhitelll994 
iBaugh. Gaztanaga. fc EfstathimJll99,^ . Both methods 
produce a particle distribution that is in unstable 
equilibrium with zero gravitational force at the location 
of each particle. In contrast, a Poisson distribution 
would immediately evolve into nonlinear structures, even 
though all three methods should recover the continuum 
case in the limit of large N. It seems better to use lattice 
or glass pre-initial conditions because they contain no 
gravitational forces, but it is not obvious that they solve 
all the problems of the Poisson pre-initial distribution 


llCIabrielli et al.ll200.3ll . 

Another method is to modify the potential of a one- 
component plasma, which will produce a particle distri¬ 
bution with the desired statistical properties of a cosmo¬ 
logical s imulation without first generating pre-initial con- 
ditions l|Gabrielli et al.l207i^l,Tovce. Levesque, fc MarcosI 
l2f)f)l . So far this method is too computationally expen¬ 
sive to be practical for large A^-body simulations. 

For simplicity I use a lattice for pre-initial condi¬ 
tions. A lattice is just an infinite series of delta func¬ 
tions at equal spacing, notated as IE, and pronounced 
“shah” (iBra.ceweHl 1200(11 . Generally, like the Gaussian, 
the in function is its own Fourier transform. Particu¬ 
larly, the Fourier transform of our dimensionless den¬ 
sity particle lattice of spacing Ai = is just V 

for k = (z,j, fc)2fcNy where {i,j,k) G all integers except 
(0, 0,0), and zero otherwise, where the “natural” Nyquist 
frequency IcNy = tt/Ax. Thus A-body simulations must 
strike a compromise between having a large enough box 
size L while ensuring that all the modes of interest are 
sub-Nyquist, so that the lattice pre-initial conditions con¬ 
tribute no power in this regime. 

Given a pre-initial particle distribution representing a 
completely homogeneous universe, the next task is to 
perturb the particles to produce a density field that re¬ 
produces the cosmologically relevant expected statistical 
properties. I will go into some detail here because it will 
provide the framework for the algorithm of El and also 
because most of the available literature on initial condi¬ 
tions is either quite abbreviated or, in some instances, 
incorrect. 

In the conventional method, the power spectrum 
of the matter fluctuations in the universe, deter- 
mine d from linear Boltzmann in tegrators such as cmb- 
fast_jl2ini^53iHiS ^1199611 or the earlier fit by 
iBardeen et alJ ill98(11 . is sampled on exactly those values 
of k which are the abscissae for the finite-volume power 
spectrum. The power spectrum is usually calculated in 
linear theory and extrapolated to the present epoch. As¬ 
suming a ACDM cosmology, in first order Eulerian per¬ 
turbation theory all modes evolve independently, so the 
power spectrum can be scaled back to the initial epoch 
for the cosmological simulation via the growth function 


D{a) = 


SBihAq 




( 2 ) 


10 


llDodelsonl200,‘1. equation 7.77) where dots denote deriva¬ 
tives with respect to proper time. The time derivative of 
the growth function is 


A(a) = 

2aa 


^ 3D 2n],D 
a 


(3) 


We will mostly use the normalized growth function 
D{a) = The Hubble parameter at a given epoch is 


- = (4) 

a 

where Hq, Bm, and are parametrizations of the cos¬ 
mology and Bk = 1 — Bm — Ba, and radiation is ignored 
here. 
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The power spectrum, extrapolated to the present 
epoch in linear perturbation theory via the growth equa¬ 
tion, is normalized to a given value of as such that 

= ^ J P{k)W'^{k)k'^dk (5) 


where W (x) is a top-hat spherical window function of ra¬ 
dius R = 8 /i“^Mpc, normalized to integrate to 1, whose 
Fourier transform is 


W{k) = W{k) 


——^ (sin fci? — kR cos kR). 
(kR)'^ 


( 6 ) 


For a Gaussian random field, a single realization should 
populate a given mode i5(k) with a value whose phase is 
random, and whose magnitude is drawn from a Rayleigh 
distribution^’^ 


PTS,{6)dS = -^e-^dS (7) 

where the variance cr^ = VP{k)/2 ensures that 

(|i5(A:)p) = VP{k). This is equivalent to drawing 
both the real and imaginary parts of (5(k) from Gaus¬ 
sian distributions with va riance a lso a^ = VP(k)/2 

S a fc Padmanabhaiil Il997t iKlvo in fc HoltzmanI 

The density field is real, so i5(k) must satisfy the usual 
Hermitian constraints. Particular attention must be paid 
to the eight modes at the corners of the S array which 
are forced to be real (those modes where each Cartesian 
component of k is either 0 or kjssy)- The value of 6 at 
these k should be drawn from a (real) Gaussian while still 
satisfying (i5(k)) = 0 and (|i5(k)p) = VP{k). This issue 
is especially important for the k = 0 mode (if P{0) ^ 0, 
as we will consider in fQ because there are no other 
modes with similar amplitudes of k to dilute any error. 

Obviously k cannot go to infinity in a computer pro¬ 
gram, so a Nyquist cutoff is introduced which can be 
of order the natural Nyquist frequency tt/Ax given by 
the particles. The field i5(k) can then be inverse Fourier 
transformed into the real-space density field, defined 
on a Eulerian grid. To convert the Eulerian density 
into a Lagrangian representation of p articles, the Zel- 
dovich approximation is widely used llZeldovichl 119701: 
lEfstathiou et al.lll98fill : 


X = q-f’®'(q) 

(8-a) 

X = 4'(q) 

(8-b) 

-V-’®' = (5(q) 

(8-c) 


^ Consider temporarily the limit L ^ oo. The spacing of modes 
becomes infinitesimally small and the actual distribution used to 
populate magnitudes doesn’t matter because the central limit the¬ 
orem ensures the Gaussian field will retain its statistical properties 
if there are many modes to average over. The distribution of mag¬ 
nitudes could simply be P{5)d5 = yJVP(k))d5 for example. 

Now as L is reduced to a finite value again, the value of ^(k) in 
an infinitesimally small bin in fc-space, consisting of an average of 
power from many modes in the L —> oo case, will be drawn from 
the Rayleigh distribution. 

^ Some authors mistakenly state that the magnitude should be 
drawn from a Gaussian distribution. 


where q is the initial position of the particle, and x is 
the hnal position. Eauation l8-cl is meant to apply at any 
epoch, and since 5 scales with time as the growth fac¬ 
tor (equation [3), ^(q) = D{a)^o{q) if is the linear 
gravitational field at the present epoch.The velocity 
equation becomes 


D{a) 

= DO) 


(9) 


with D(a) given by equation 0 These are comoving ve¬ 
locities; multiply by a to get the normal physical veloci¬ 
ties.® 

In practice, it is convenient to use Fourier techniques 
to perform differentiation: e.g., d/dxi —> iki, ^ 
—The usual technique is to use Poisson’s equa¬ 
tion in Fourier space to convert i5 to a potential, and 
then to convert the real-space potent ial to the gravi¬ 
tational field ^ via finite-differenci ng (lEMathknLet_aU 
1 198, ’ll: iBaala fc PadmanabhanI n~997l iPower et al.l 1200, “lir 
However, if one has sufficient computer memory, I see no 
reason not to convert S to the gravitational field directly 
via 

’^(k) = p{k) (10) 

because the density field in the box is a continuous field 
with discrete particles, in contrast to a discrete field such 
as a Eulerian grid. Similarly, the disc rete forms of the 
Green’s functions are sometimes u sed llEfstathiou et alJ 
1198111 l^gla fc Padmanabhanlll997!l . but again, they are 
not necessarily relevant to Lagrangian cosmological sim¬ 
ulations. 


3. THE FAILURE OF THE CONVENTIONAL METHOD 

Sampling the power spectrum (or the field i5(k)) leads 
to a gross underestimate of the real space variance (e.g., 
as), the correlation function, and nonlinear effects. It is 
well known that the two-point correlation function and 
power spectrum are Fourier transform pairs. But sam¬ 
pling the power spectrum is equivalent to multiplying it 
by a in function of spacing Ak = 211 jL. By the convo¬ 
lution theorem, this sampling is equivalent to convolving 
the continuous correlation function with a III of spacing L; 
this is a simply a manifestation of aliasing, e xcept in rea l 
space instead of the more usual fc-space. As Ip^ iIT^ 
pointed out, the statistical properties of real space, such 
as the correlation function and mass variance in spheres, 
can be grossly misrepres ented. _ 

I have used cmbfast (|Seliak h Zaldarriagal 1199^ to 
generate a power spectrum at the present epoch with 
the cosmological parameters llm = 0.27, IIa = 0.73, = 

ilm —f^CDM = 0.044, Ho = 71 kms“^ Mpc”^ and spectral 
index Us = 0.97, normalized to ug = 0.84. 

^ A simple check of the Zeldovich appoximation is offered by the 
approximation <5(k) = -E ^ exp (—-ik ■ [q + ^(q)]) 

^'(q)) exp (—?k ■ q) so that (5(k) = —ik - 4^(k) for small k - 4^. This 
is just the Fourier transform of emiation 18-cl 

Most authors neglect to mention that the growth function 
should be normalized to its present-day value, an acceptable over¬ 
sight only when D{1) = 1, as in Einstein-de Sitter cosmologies. 

^ Also mind the convention that cosmological distances are usu¬ 
ally expressed in units of /i~^Mpc but velocity units usually do 
not contain the factor of h. 
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Fig. 1.— Aliased correlation functions for L = 
(100,200,400,800) /i“^Mpc (thin lines). These are obtained by 
sampling the power spectrum on a grid, inverse Fourier transform¬ 
ing, and collapsing angular dependence. Because the 3-D correla¬ 
tion function array is cubical, information about ^(r) can be ob¬ 
tained for r G (T/2, -^L) from the “corners” of the ^(x) array. 
This information about ^(r) can be thought of as arising from cor¬ 
relations between the density field in opposite corners of the <5(x) 
box. The grid sizes were 128^,128^,256^, and 512^ respectively, 
which are large enough that this plot would look little different 
in the limit of infinite grid size. The thick line is the continuous 
inverse Fourier transform of the cmbfast power spectrum. 


Figure [n shows correlation functions obtained by sam¬ 
pling this power spectrum with spacing AA:, perform¬ 
ing an inverse fast Fourier transform to get a discrete 
3-dimensional correlation function, and then averaging 
values of ^(r) in the grid for each unique r. These 
aliased correlation functions are representative of the 
correlation function one would measure in a density 
field after generating initial conditions with the conven¬ 
tional method of (neglecting linear growth). Also 
shown is the infinite-volume correlation function, ob¬ 
tained by an exact Fourier transform of the power spec¬ 
trum (see equation 113; the baryon bump is visible at 
r « 120 /i“^Mpc. It is easy to see that box sizes as small 
as 50 or 100 /i“^Mpc have profound effects on the corre¬ 
lation function at virtually all length scales. Note that 
since F’(O) = 0, the integral of the correlation function 
over the box must equal zero; this explains the fact that 
the zero-crossing tq of the correlation function occurs at 
substantially smaller r than in the continuous case, es¬ 
pecially for box sizes L < 2ro. 

This problem is frequently stated as the finite 
value of L cutting off pertur bations with wave¬ 
lengt hs g reater than L fcf. FO elb fc Bertschine^ 
19^ _ iBaugh. Gaztaiiaga. fc EfstathioiJ Il995t 

T 2 rmen_jj_jIgr^g(,dfinggH Il99fit iKnebe fc Dominguez! 
200.‘It iBagla fc R.avll2004D . Actually, the finite value of 
L cuts out all modes k that do not lie on the lattice 
n of spacing Afc. However, the effect of sampling may 
be more pronounced at super-L modes because the 
three-dimensional power spectrum has a cusp at k = 0 
(i.e., F’(k) cx for small k where n « 1 for currently 
favored near scale-invariant cosmologies). 

A preliminary understanding of the effect of super-A 
modes on the mass variance in spheres can be obtained 


by integrating 

4 = ^ P{k)\W{k)\^edk (11) 

for vari ous values of A, whe r e Wi k) is given by equa¬ 
tion El iGelb fc Bertsching^ lll994ll found that the un¬ 
derestimate of (Tg for A = 50 /i“^Mpc was about 10% for 
an Einstein-de Sitter cosmology. Since currently favored 
ACDM cosmologies have more power on larger scales (for 
a given power spectrum normalization), this underesti¬ 
mate would be even more severe if calculated with the 
cmbfast power spectrum used in this paper. 

However, the integral of equation is not quite the 
correct method for determining the effect of super-A 
modes. Since the real-space box is a periodic finite vol¬ 
ume, the Fourier convention of equation^should be used. 
Then the equation for o-g is accurately modified to 

^1 = ^ E Ak)imk)P (12) 

kG(i,j,fe) Afc 

and the contribution to (t| from super-A modes is given 
by P(0)/E, scaled by A)^(a) if earlier epochs are consid¬ 
ered. I will return to this point in a but first I must 
discuss how to get the value of the power spectrum at 
k = 0 by convolving the power spectrum. 

4. CONVOLVING THE POWER SPECTRUM 

Fr om the point of view of the lPress fc Schechted lll974ll 
and iHamilton etalJ (11991|) theories of nonlinear struc¬ 
ture formation, iPenI lll997ll argued that the real space 
statistical properties should be accurately modeled inside 
the simulation box. This is conceptually done by mul¬ 
tiplying the real-space correlation function by the win¬ 
dow function representing the box geometry; i.e. setting 
^{xi) = 0 for l^il > A/2 where Xi is a Cartesian coordi¬ 
nate. From the convolution theorem, this is equivalent to 
convolving the power spectrum with the Fourier trans¬ 
form of the window function. Finally, the real-space box 
is made periodic by convolving it with a 11 of spacing 
A, which is equivalent to sampling the power spectrum 
with a in of spacing Afc. It is the crucial eonvolution 
of the power spectrum before sampling that ensures that 
the real space statistical properties inside the box will be 
correct. The mass variance in spheres of radius R will be 
accurately modeled for box sizes A > AR. 

I found that it was easiest to convolve the power spec¬ 
trum by first inverse Fourier transforming the continuous 
power spectrum into the correlation function, truncating 
the correlation function at A/2, and then Fourier trans¬ 
forming back. The relevant isotropized Fourier transform 
equations are: 

^ ^ Jo 

sin kr 

Phik) =Att f{r) —-—(13-b) 

Jo 

Upon sampling this convolved power spectrum P^fk), the 
correlation function of the resulting real-space field will 
be exact in a sphere of radius A/2, and the mass variance 
of spheres will be exact for radii i? < A/4. It is concep¬ 
tually possible to ensure that the correlation function is 
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exact in the entire cubical box instead of just its inscribed 
sphere, but this would be numerically more challenging, 
and it is unclear that it would be of any benefit because 
of the anisotropy of the box on large scales. 

The sampling of the convolved power spectrum at a 
given k can be understood as taking a suitable weighted 
average of the power in a bin of characteristic size Ak^ 
centered on k liPenlll 1)9711 . This will have a particularly 
pronounced effect on the k = 0 mode, due to the cuspy 
nature of the power spectrum at k = 0. Another way 
of interpreting the sampling of a convolved power spec¬ 
trum is that it is the Fourier equivalent of convolving the 
truncated correlation function with a IH of spacing L, i.e., 
making it periodic (and introducing aliasing error in real 
space if the correlation function was not truncated!). 

Convolved power spectra for various values of L are 
presented in figure |21 The original power spectrum is 
again given by cmbfast. The resulting correlation func¬ 
tions for these power spectra are shown in figureOl which 
can be compared directly to figure Q] As expected the 
correlation functions are exact within L/2. 

5. THE DC MODE 

In order to get the real-space statistics of the density 
field correct within the box, it is clear that the power 
spectrum must be convolved as described in m before 
sampling. The generation of initial conditions can then 
proceed as usual, as in fj21 with one additional nontriv¬ 
ial correction: Pl( 0) is no longer assumed to be zero. 
Figure 01 shows what this power is for various box sizes; 
ideally one would run a simulation with L oo and 
Pl{0) —*■ 0.® The relevant Green’s function (equation^)) 
becomes undefined at k = 0. What should one do about 
this DC mode? 

Note that the amplitude of the DC mode is drawn from 
a Gaussian, so it will be necessary to run an ensemble of 
simulations with different random number seeds. 

But first, consider a single realization. The DC mode 
corresponds to the imposition of a constant overdensity 
at present epoch Aq drawn from a Gaussian distribu¬ 
tion of variance Pl{0)/V. Since it is a large-scale mode, 
its amplitude will evolve according to the linear growth 
function (equation^ : A = DAq. For small A, one could 
evolve an A-body simulation ignoring the DC mode, 
and then add the contribution from A to an evolved 
simulation by violating mass conservation and scaling 
the mass of each particle by 1 -|- A. Alternatively, one 
could rescale the box of an evolved simulation, effectively 
p ushing the particles as in the m o de-adding p rocedure 
of iTormen fc Bertschinge j ijlQQlij) . iCold ijlQQTf) showed 
that, in addition, one should rescale cosmic time in or¬ 
der to more accurately emulate long-wavelength mode 
coupling. This coupling of long-wavelength modes (and 
by extension, the DC mode) to all other modes can be 
thought of as follows. A large-scale overdensity effec¬ 
tively modifies the local value of Dm (and Da) in its do¬ 
main of the universe. That is, structures in overdense 
regions grow faster. I discuss this quantitatively below 
and in Appendix IXI 

It would be even more accurate to incorporate the 
effects of DC mode coupling into the simulation it- 

® I am assuming Ne wtonian gravity eve n for the largest boxes, 
as is usually done l'e.B'..lEvrard et al.l l200j) . 


self. Since A > 0 corresponds to a box which ex¬ 
pands slower than the A = 0 case, and vice versa for 
A < 0, we just need to determine abox(0j where t is 
cosmic time, so that the modified equation of motion, 
X -I- 2xabox/abox = ~'^4>/<^hox’ accurately encapsulates 
the effects of DC mode coupling. I will use the subscripts 
“box” and “uni” to refer to any parameter that relates to 
a simulation with Aq ^ 0, and Aq = 0, respectively (i.e., 
Ouni is the scale factor of the universe, Dm, uni is given by 
the cosmological model but Dm, box may differ, etc.). 

Since auni(t) is known from equation 0] it is suffi¬ 
cient at first to derive the relationship between Obox 
and Ouni- From a Eulerian viewpoint, writing the ra¬ 
tio of the densities of the box and the universe results 
in Obox/ouni = (1 + A)~^/^. From a Lagrangian view¬ 
point, Obox/fluni = (A/2 — 'I')/(L/2) where T is a small 
displacement equal to the force field by the Zeldovich 
approximation: V • SF = A. Considering a test particle 
at coordinates (A/2,0, 0) where the origin coincides with 
the center of the box, T = LA/6. Therefore, 

^ = l-iA. (14) 

Uuni 3 

I will use this relationship between Obox and Ouni on the 
grounds that Lagrangian perturbation theory is, gener¬ 
ally sp eaking, more accurat e than Eulerian perturbation 
theory llBouchet et alJl99,^ . Obviously both approaches 
agree as A —> 0. 

Now abox(i) is known, so the A-body code could be 
modified to reflect this change in the equation of motion. 
However, one could avoid rewriting A-body code (and 
gain insight into the physics of the DC mode) by modi- 
fiying the cosmological parameters of a given realization 
with A 0 as follows: 


-^O.box — -^O.uni ^ , , 

(15-a) 

^m,box “ ^m,uni(l + 0) 

(15-b) 

^A,box — ^A,uni(l + 0) 

(15-c) 

D 


(16) 


I provide the derivation in Appendix m This is the 
prescription for running an A-body simulation with a 
nonzero DC mode. Modify the cosmological parameters 
as appropriate for a given realization, run the simula¬ 
tion as normal, but map the cosmic time according to 
equation 01 That is, if you are interested in the simu¬ 
lation at epoch ^uni, run the simulation until Zbox given 

by 1 -k 2box = (1 + -2uni)/(I - A) An/3). _ 

In their numerical simulations. iFrenk et alJ llI988t) used 
essentially the same scheme as outlined here to treat 
the DC mode. They ran sets of simulations with the 
DC mode equal to {0, ±1} times the rms overdensity, 
i.e. \/A’l( 0)/E, instead of drawing the DC mode from a 
Gaussian as in this paper. They also synchronized their 
ensemble at z = 6 with every realization having the same 
box size (2 Mpc), thus varying the integrated mass in the 
box. In contrast, the present approach synchronizes the 
box size at 2 —> oo (see Appendix 0 and uses the same 
box-integrated mass among the realizations. 
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Fig. 2.— The cmbfast power spectrum (black line) and convolved power spectra for L = (50,100,400,800,1600,3000) h ^Mpc (colored 
lines). The bottom panel shows the convolved power spectra normalized by the cmbfast power spectrum on a linear scale. (The 1600 and 
3000 /i~^Mpc cases are not plotted in the bottom panel.) The ringing in the convolved power spectra has wavelength 47r/L. Note that 
each convolved power spectrum is intended to be sampled with a specific interval Ak = 277jL (in three dimensions), but I have plotted 
them as continuous functions for illustration. 



Fig. 3. — Same as figure^ except that the appropriate convolved 
power spectra are sampled for L = (100,200,400,800) h“^Mpc. 
As intended, the correlation function matches the L oo limit for 
r < L/2, and is 0 for r > L/2. 



Fig. 4. — Power in the DC mode vs. box size L (top panel) and its 
dimensionless version, the characteristic dimensionless overdensity 
corresponding to Pl( 0) (bottom panel). I used the cmbfast power 
spectrum normalized to as = 0.84. 


6. ENSEMBLES OF SIMULATIONS 

To test the new method for generating initial condi¬ 
tions, I have run several ensembles of simulations. I will 


refer to initial conditions generated in the traditional 
way, by sampling the power spectrum, as P-sampled 
IC’s. Initial conditions generated with the new method 
of convolving the power spectrum before sampling can 


































be called ^-sampled IC’s because the correlation func¬ 
tion will be accurately reproduced in the box7 

An ensemble of simulations is a set of A^-body simula¬ 
tions with the same cosmological parameters and techni¬ 
cal parameters (e.g., N,V), but with different random 
number seeds. For simulations running on ^-sampled 
IC’s, the effective cosmological parameters of a given re¬ 
alization within the ensemble change as described in 
but the realization still represents a part of the universe 
that the ensemble is intended to model. 

For each ensemble, I generated initial conditions to 
second order in Lagrangian perturbati on theory (2LPT, 
lScoccimarrdll998tlBmicliet et al.iri995ll for the following 
output redshifts: 

log [1/(1 + zuni)] = -3.0, -2.8, -2.6, 

-2.4,-2.2,-2.0,-1.8,-1.6. (17) 

I used the last of these outputs as initial conditions to the 
A-bo dy code, tree particle mesh ITPM. rRode fc Ostrikerl 
unni. The output redshifts from TPM were 

1/(1+ ^u„i) = 10-^■^ 10-1-2, 

0.10,0.15,0.20,0.25,-••1.00 (18) 

All the ensembles presented herein comprised 100 re¬ 
alizations of the cosmology flm = 0.27, flA = 0.73, h = 
0.71. The only parameters varied in this work are L, A, 
and whether I used P- or ^-sampled IC’s. I used the same 
100 random number seeds for all ensembles, and because 
of the way the Fourier modes are populated in the initial 
conditions, the same large-scale structures, scaling as L, 
develop in all ensembles for a given random seed. 

Note a subtle change in terminology at this point. In¬ 
stead of inputting parameters like as,P{k) into initial 
conditions, I will now present and analyze measured val¬ 
ues for (78 j P{k) E^nd ^(r) from the particle distributions 
resulting from the combination of 2LPT and TPM. 

7. THE MASS VARIANCE IN 8 ft-iMpc SPHERES cts 
7.1. Method 

To calculate erg for a given epoch of a single realiza¬ 
tion, it is possible to use equation d if the power spec¬ 
trum is known. The discrete Fourier transform, for which 
the well-known fast Fourier transform is an implementa¬ 
tion, introduces aliasing error into the power spectrum. 
Even if the continuous Fourier transform is used (equa¬ 
tions^, it is impossible to have fc —> oo in a computer 
program. It seems better to be able to control the er¬ 
rors easily by using a real-space Monte Carlo integration 
scheme: 


^8 = ^E(<5(xr)*VF)" (19) 

<5(xr)*IT = (20) 

^sphere 

where Ar random points Xr within the box volume are 
chosen to serve as integration abscissae. I4ni is the vol¬ 
ume of the box if Aq had been zero. The window function 
IT is a spherical tophat of comoving radius 8 h“^Mpc 

^ The correlation function is not actually “sampled” in the usual 
sense. 



(8 Mpe/h) / L 


Fig. 5.— The value of erg lattice a® a function of Rj/\x, where 
the mean interparticle separation A.x = L/N^t^ and the radius of 
the spherical tophat is i? = 8 /i~^Mpc. Values are accurate to 1%. 

and volume Tpherel the number of particles within that 
sphere (at location x^) is simply Agphere- The error on 
(tI can then be estimated as 


= 


'((^*IT)4) - {{5*wy 


At 


( 21 ) 


llPress et al lIUll §7.6) where brackets indicate an av¬ 
erage over At random points. Unless otherwise noted, 
all values of ct| in this paper have been calculated to 2% 
accuracy (for individual realizations), which corresponds 
to approximately 1% accuracy on erg itself. 

Because (t| is a volume average (assuming ergodicity), 
the value of o-| for the ensemble is a weighted average of 
the values of cTg from each realization, and the weight of 
each realization is 

Let CTg, lattice be the value of erg from an unperturbed 
lattice of particles, and let crg,poisson be the expected 
value of CTg from a Poisson distribution of particles. T hese 
two quantities will be important to the discussion in o 
below. 

It is easiest to calculate CTg lattice with a Monte Carlo 
integrator, using equation Figure |3 shows the re¬ 
sult, which only depends on the ratio of the radius of the 
spherical tophat R = 8 h”^Mpc to the spacing between 
particles Ate = L/N^^^. Naturally, as Ax decreases, the 
effects of particle discreteness become less severe, and 
CTg,lattice decreases. The wiggles in figureElcan be thought 
of as interactions between the spherical symmetry of a 
sphere and the translational symmetry of a lattice as the 
relative scales change. 

A Poisson distribution of particles results from plac¬ 
ing particles randomly throughout the box. It is well 
known that a Poisson distribution has F’(k) = V/N for 
all k. Note that V/N is the expectation value, and power 
spectra of realizations and finite ensembles will fluctu¬ 
ate about this value. The expectation value of Ug poisson 
should properly be calculated from equation El because 
the real-space box is a periodic finite volume. However, 
we can avoid evaluating the infinite sum by using the fol¬ 
lowing trick, based on calculating CTg in real space. Imag¬ 
ine a box of size L > 4i? into which we will throw down 
particles, randomly placed. As we throw down particles. 
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imagine convolving the particles with the R = 8 /i“^Mpc 
spherical tophat so that each particle becomes a diffuse 
sphere. Since the box is periodic we can assume the first 
particle lands in the exact center of the box without loss 
of generality. The second particle may land so that the 
distance between the two particles is either less than 2R, 
or greater than 2R. In the former case, the “domains 
of influence” of the two spheres overlap and it is exactly 
this “interaction” that our real space erg calculator must 
identify. In the latter case, the distance between the two 
particles is irrelevant. In other words, L could have been 
twice as large and the distance between the two parti¬ 
cles could have been twice as large, and our ug calcula¬ 
tor would not have been able to tell the difference. (The 
introduction of more particles just generalizes to a super¬ 
position of many two-particle “interactions.”) Therefore, 
as long as L > 4i?, the size of the box does not matter 
for our (Tg calculator. If L < AR, then the second particle 
could “interact” with the first particle on both sides due 
to periodicity, and this double-interaction would not be 
a Poisson process. The above argument implies that the 
infinite sum of equation d must yield the same (exact) 
answer as the continuous integral 

1 r°° 

(22) 

for a Poisson distribution in a box of size L > AR, where 
W{k) is given by equationEl The value of the integral is 
3V/AttR^N, so cTg,Poisson = 0.1193 for L = 100 /i“^Mpc, 
N = 32^. Strictly, for a Poisson distribution the number 
of particles in the box N fluctuates on the order of VN, 
but I will use the term “Poisson” even when N is fixed. 
The effect of using exactly N particles is to set the DC 
mode to zero, and the missing contribution to (t| from 
this single mode is 1/iV, too small to fret about in this 
work. 

7.2. Results 

Figure El shows erg for two ensembles with L = 
50 /i“^Mpc, N = 32^ differing in whether they use P- 
or ^-sampled IC’s. Since a plot with all 100 realizations 
would be too confusing, I have plotted the range of ug 
for the 100 realizations in light green. The 15th and 84th 
percentiles are bounded by the purple region, which ap¬ 
proximately indicates the 68.3% confidence region for in¬ 
dividual realizations. I used a bootstrap method to deter¬ 
mine the 68.3% confidence region on the actual ensemble 
average, plotted as the black region (more optimistically, 
the thick line). The ensemble average can be compared 
to the theoretical expectation Ug^theory = D{a) x 0.84, 
plotted as the lower red line. It is evident that a sys¬ 
tematic bias exists in the value of erg for the ensemble in 
the P-sampled case; it is too low. The ^-sampled case 
corrects this bias as expected. I will discuss this effect 
more quantitatively below. 

The behavior of the ensemble value of erg at early times 
is well approximated by (er^^ti^g -k cr8,theory)^^^ plot¬ 
ted as the middle red line in figure El Why does the 
contribution from the lattice add in quadrature? Recall 
that the power spectrum of the lattice of particles is zero 
everywhere except for modes with each component an 
integer multiple of 2fcNy (excluding k = 0). Applying 
a small displacement field to the particles results in a 


power spectrum for the particles that accurately mod¬ 
els the input power spectrum for sub-Nyquist wavenum¬ 
bers. If the displacement field is small, then the lattice 
part of the power spectrum will be barely modified for 
k ~ 2A:Ny, but will be more greatly modified for very 
large k. From equation 1121 the measured value of cr| 
can be decomposed into two sums. The first sum, over 
all sub-Nyquist wavenumbers, measures Ug t^gory from 
the information in the displacement process. The sec¬ 
ond sum, over all super-Nyquist wavenumbers, approxi¬ 
mately measures cr| lattice because the weight of |lF(k)p 
is greatest for those modes which are least affected by 
the displacement, i.e., the modes k ~ 2fcNy. 

The upper red line in figure El shows ((t| Poisson + 
trl theory) fr i® dear that Poisson noise is not a good 
model for the underlying discretization of matter into 
mass particles for most of the logarithmic range of the 
simulation. In fact, vestiges of the initial lattice can be 
seen in these simulations to z ~ 2 (and even later times 
for simulations with larger Ax), so there is no reason 
to model the discretization of matter into particles as a 
Poisson distribution before these epochs. 

Figure [3 shows the evolution of erg for both P-sampled 
and .^-sampled IC’s for five different combinations of 
technical parameters {L,N). In this figure all curves 
are normalized to (ul^theory + so to the ex¬ 

tent that the simulations evolve according to linear the¬ 
ory the black region should match the lower red line 
0-8, theory /(cTg, theory + o"!, lattice) • The middle two pan¬ 
els, representing the case P = 50 /i“^Mpc, N = 32^, 
show the same information as figure El for comparison. 
The two rows above it have N = 32^ as well, but 
L = (100, 250) /i“^Mpc so that ug,lattice is large. The 
bottom two rows of the figure show the effect of increas¬ 
ing L while keeping Ax fixed. As expected, in the P- 
sampled case, the discrepancy between ug and erg, theory 
(for intermediate times) is reduced for larger P, and com¬ 
pletely eliminated in the ^-sampled case. It also appears 
that nonlinear evolution of ug at late times is enhanced 
using the larger box size. This is presumably due to the 
greater density of modes available for coupling, and can¬ 
not be solved using ^-sampled IC’s with small box sizes. 
However, it is fair to claim that using ^-sampled IC’s 
solves the first order problem of the discrepancy between 
(Tg of the ensemble and erg, theory 

Let’s now take a step back to examine the effect of 
P-sampled vs. ^-sampled IC’s on individual realizations. 
FigureElshows the values of erg for the 100 realizations of 
P = 50 /i“^Mpc, N = 32^ at z = 38.81 (left panels) and 
z = 5.67 (right panels). For this figure I just used the 
2LPT initial conditions code for both of these redshifts, 
and I calculated values of er| to 0.5% accuracy. In the 
upper panels, note that the values of erg for the ^-sampled 
case (blue filled circles) tend to be greater than those for 
the P-sampled case (green open circles), and that the 
discrepancy is greater for larger |Ao| (of the .^-sampled 
IC’s). I have also generated particle positions using the 
convolved power spectrum but artificially setting the DC 
mode to zero; the values of erg for these configurations 
are plotted as red triangles. By comparing these three 
types of IC’s, it is apparent that most of the discrepancy 
between the P- and ^-sampled IC’s comes from the DC 
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Fig. 6.— erg as a function of scale factor Uuni for -P-sampled IC’s (left panels) and ^-sampled IC’s (right panels). The ensembles had 
L = 50 h”^Mpc, N = 32^. For a given auni? I calculated erg with a Monte Carlo integrator for each realization. The light green region 
shows the range of erg for the 100 realizations. The purple region shows the 68.3% confidence region on individual realizations, and the 
black “line” shows the 68.3% confidence region for the ensemble average. The linear theory expectation of erg’s evolution is plotted as the 
bottom red line. The middle red line is the theoretical value with the expected lattice contribution added in quadrature; the top red line 
is the theoretical value with the appropriate Poisson contribution added in quadrature. The bottom panels plot the same information as 
the top panels, except using a linear scale, which makes clearer the discrepancy between the theoretical and calculated evolution of erg in 
the P-sampled case, and the agreement in the ^-sampled case. In the bottom panels I have also plotted, with blue lines, ten of the 100 
realizations to show qualitatively individual errors due to the Monte Carlo integrator. I chose the same ten random seeds for both the 
P- and ^-sampled cases. These individual evolution tracks also show qualitatively that a realization with a small value of erg early in the 
simulation has a small value later in the simulation as well. 


mode itself. 

The lower panels in figure |H1 show the difference in 
cr| between the ^-sampled and P-sampled IC’s, = 

f^l.C-sampied “ '^l.p-sampied' ^s I claimed above, the 
dominant difference between and P-sampled IC’s man¬ 
ifests itself in the DC mode, (t| is approximately the 
contribution to ct| from the DC mode. From eauation ll2l 
the expected contribution from the DC mode is just 
= P^(a)PL(0)/V = (P(a)Ao)^, plotted as a line. 
The expected contribution agrees well with tT| , espe¬ 
cially for higher redshift when, presumably, mode cou¬ 
pling has less effect. 

Now we can answer the question: what is the effect of 
using ^-sampled IC’s instead of P-sampled IC’s on the 
ensemble value of erg? For early times, the weighting 
of each realization (by properly required) is not 

important, so the ensemble contribution cr| ^ is an 


integral over the realization values cr| 


8,DC' 


CTi 


8.DC,e - 

I 




<^I.Dce (23) 


where Pl{0)/V is the variance of the DC mode given by 
the convolved power spectrum, and Aq is the dimension¬ 
less DC mode of_a single realization (at present epoch). 
Since cr| r; (P(a)Ao)^, the expected contribution to 
the ensemble average CTgpQg = Pl{0)/V at the present 
epoch. The error in ug due to the DC mode when using 
P-sampled IC’s is approximately /V)/(^1 

where Pl( 0)/V can be read off of figure 0 Figure 0] 
uses a normalization of erg^theory = 0.84 (which can¬ 
cels out in the estimate of this error). The error is 
about 6% for P = 50 h“^Mpc, and reduces to 1% 
for L = 120 ft,“^Mpc. This is not to say that a P- 
sampled simulation of the latter size is accurate to 1% 
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Evolution of as of ten ensembles, normalized to the linear theory prediction plus the expected lattice contribution ((t| - 

The light green, purple, and black regions, and the three red lines are as in figure El 


Fig. 7. 

0-2 fi/a 

^8,lattice/ 

in all respects; note that the characteristic dimension¬ 
less overdensity for L = 120 /i“^Mpc is about 12% (foi 
L = 50 /i-^Mpc, 30%). 

8. THE CORRELATION FUNCTION ^{R) 

8.1. Method 


The well-known formula for the correlation function is 

e(x) = (<5(x')<5(x-l-x')) (24) 

where the brackets indicate an ensemble average, or space 
average assuming ergodicity. But first, we will need the 
correlation function of a single realization. Substitute 
6 = p/pq — I, where po is the ensemble (universe) average, 
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Fig. 8.— The upper panels show 100 individual realizations’ values of crs for L = 50 /i“^Mpc, N = 32^ for z = 38.81 (left panels) and 
^ = 5.67 (right panels) with particle positions generated with the second order Lagrangian perturbation theory initial conditions code. The 
green open circles are for P-sampled IC’s; the blue filled circles are for ^-sampled IC’s, and the red triangles are for ^-sampled IC’s with 
the DC mode artificially set to zero. All three cases used the same 100 random number seeds. A thin black line connects the range of 
values for each realization. The values of Aq from the ^-sampled case serve as the abscissae. Values of erg are accurate to 0.25%. The lower 
panels plot the 100 realizations’ values of cr| qq = a-g j_sajupled “ ^8 P— sampled' expected value of = (D(a)Ao)^ is also shown. 


and isotropize to obtain the formula for the correlation 
function of particles in a box of size ibox^ 

^uni ^uni (A^shell) 
t^box / h^hell 

-2(—] +1. (25) 

\ flbox / 

Here A^jbeii is the number of particles in a thin spherical 
shell of radius r and volume 14hen surrounding a single 
particle, and (iVsheii) is an average over all the particles 
in the box (or a random subset). 

For L/2 <r < the “shell” used in equation EKl 

becomes the intersection of a spherical shell -with the box. 
In this way, it is possible to obtain information about the 
correlation function from particles on opposite corners of 
the box, but it must be stressed that these are the scales 
on which the box is anisotropic so it is unlikely that the 
correlation function on these scales should ever be used 
in rigorous analysis. 

The correlation function of the ensemble is a simple 
weighted average of the correlation functions of the indi¬ 
vidual realizations, and the weight of each realization is 


“box- 

8.2. Results 

Ensemble correlation functions are presented in fig¬ 
ure El for several output redshifts for two ensembles, dif¬ 
fering in whether they used P-sampled or ^-sampled IC’s. 
These ensembles had L = 100 Ii“^Mpc and N = 64^. I 
binned the correlation function into 0.3 /i“^Mpc bins. 
Note the large systematic offset between the measured 
correlation function and linear theory prediction in the 
P-sampled case, which is largely corrected in the ^- 
sampled case for r < P/2, as intended. The “ringing” in 
the correlation function at early times is due to the initial 
lattice and is erased cleanly as the simulation evolves. 

9. THE POWER SPECTRUM P{K) 

9.1. Method 

To eliminate aliasing error, I used the brute force 
Fourier transform (eauation ll-aH to determine the power 
spectrum of each realization. I chose the following ab¬ 
scissae (k = {i,j,k)Ak): 

(uj,fc)G±{0,l,2,3,'"16,18,20,'" 

32,36,40" •64,72,80,"-128} 


(26) 


























13 


o 

CL 


60 

40 

20 

0 

-20 

-40 

-60 

15 

10 

5 

0 

-5 

-10 

-15 

2 

1 

0 


-2 


P —sampled ^ — sampled 



0 20 40 60 80 20 40 60 80 100 

r [Mpc/h] r [Mpc/h] 

Fig. 9.— Correlation functions for the L = 100 /i“^Mpc, N = 64® ensembles with P-sampled IC’s (left panels) and ^-sampled IC’s (right 
panels) for three output Uuni? labeled. The black region is the 68.3% confidence region on the ensemble correlation function; the purple 
region is the approximate 68.3% confidence region on individual realization correlation functions, and the light green region bounds all 100 
realizations’ correlation functions. In the upper four panels, five individual correlation functions for example are overplotted in blue (the 
random number seeds chosen match across each panel). The red line is the linear theory prediction. As in figures [Tl and information 
about the correlation function is available for r > (nbox/^iuni)^/2 from the corners of the simulation cube, but it is unlikely that this 
information should ever be used in rigorous analysis. 
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to achieve a balance between good coverage in fc-space 
and computational feasibility. It is important to consider 
biases introduced by this sampling of abscissae. For ex¬ 
ample, the power spectrum of a lattice of 64^ particles 
has a large spike at k = (64, 64,64)AA: and is zero in 
the bins surrounding that k, so the power in that mode 
should not be taken to be representative of the power in 
the surrounding bins. 

For .^-sampled IC’s, the effective box size of each real¬ 
ization is modified, and therefore Ak is different among 
the realizations in an ensemble. The power spectrum of 
the ensemble, like the correlation function and (t|, can 
be determined from the realizations’ power spectra by 
weighting each by before averaging. (This can be 
seen from the linearity of the Fourier transform of the 
correlation function.) For the ^-sampled case, the ab¬ 
scissae of each realization shift, so I binned the power in 
fc-space before averaging. I used a binning correspond¬ 
ing to the spacing of fc in equation[2El To be consistent I 
binned the P-sampled case correspondingly. It was nec¬ 
essary to keep careful track of the number of modes at 
each fc = |k| to do the averaging correctly. 

9.2. Results 

Figure cni shows power spectra for the same two en¬ 
sembles discussed above in ^8.21 for various epochs. The 
DC mode is omitted from these plots. Not much is differ¬ 
ent between the P-sampled and ^-sampled cases, except 
for the oscillations in the .^-sampled case for early epochs 
which are subsequently washed out, and the much greater 
dispersion in normalization for the ^-sampled case at late 
epochs. The latter effect is a consequence of the mode 
coupling of the DC mode to all other modes; a box with 
a large-scale overdensity effectively emulates a universe 
with a larger Dm, etc. However, as we have seen directly 
in figure inland indirectly in figure HI the DC mode is the 
most important effect of using ^-sampled IC’s. Therefore 
the difference in ensemble power spectra (besides the DC 
mode) between P- and ^-sampled IC’s is not as dramatic 
as the difference in correlation functions seen in figure [51 

10. SUBTRACTING POISSON NOISE 

When measuring power spectra of N-body simula¬ 
tions, the Poisson noise contribution is sometimes sub¬ 
tracted to arrive at a more “true” power spectrum that 
represents the underlying con tinuous density field (e.g., 
iSmith et alJ 1200,11 : l.linel 120041) . The rationale is that 
the discreteness of particles introduces a delta func¬ 
tion in the correlation function, or that a Poisson dis¬ 
tribution is representative of a field with “no” cluster¬ 
ing. However, 7V-body simulations start with a pre¬ 
initial lattice (or glass) which is not a Poisson distri¬ 
bution of particles. As shown in figures 0 and [T] the 
Poisson contribution to erg i s a bad model until at leas t 
rat her late epochs. Similarl y. IBaugh fc Efsta thimJ (119941) 

and iBaugh. Caztanaga. fc EfstathiouTIT^iTl) pointed out 

that the Poisson contribution cannot be subtracted until 
the power spectrum evolves above the level V/N. Look¬ 
ing at the (Tg results of figures 0 and 0 it is unclear 
whether or not the Poisson contribution is a bad model 
at very late epochs, though, since the lattice is mostly 
erased by that time. 

Another insight is provided by the power spectra of fig¬ 
ure cni Note the spikes in the power spectrum at early 


epochs. These are caused by the lattice of particles and 
initially have power H at k = (L 7 , fc)2 7r/Aa:. k 0. 
Thes e spikes can be called Bragg peaks ilGabrielli et alJ 
120051) by analogy with x-ray scattering of crystals. Par¬ 
ticle evolution smooths out these Bragg peaks, so that 
at a ~ 0.2 in figure El (before the imposed perturba¬ 
tion spectrum has significant effect on high fc) the power 
becomes Poissonian at high fc. That is, the power has re¬ 
distributed to an expectation value V/N over the N bins 
in fc-space surrounding the initial Bragg peak. The inte¬ 
grated power has been neither magnified nor suppressed 
because the power at these high fc is due to discrete¬ 
ness, not evolution of a collisionless fluid. The crucial 
point is that there is no initial power in the lattice at 
k = 0, so at sub-Nyquist fc, the resulting Poisson power 
spectrum contribution does not exist, and should not be 
subtracted. The lack of the Poisson power spectrum at 
sub-Nyquist scales can be seen from the middle panels 
of figure El where there is a dip in the power spectrum 
below the Poisson level at fc ~ 2.5 fcMpc”^. However, it 
remains unclear if the discreteness power at high fc can 
migrate inwards to affect all fc by the end of the sim¬ 
ulation, and therefore if some reasonable model of the 
discreteness contribution to the power spectrum should 
be subtracted. The problem of the spreading of initial 
discreteness power to the wavenumbers of interest in the 
simulation will be left as an open question. The issue 
may be made moot merely by using negligible values of 
V/N in simulations of the future. 

11. SUMMARY 

I have desc ribed and i mplemented ^-sampled IC’s, as 
suggested by iPenI (TT^ . The advantage of ^-sampled 
IC’s is that they accurately reproduce the real-space sta¬ 
tistical properties of the universe within the box. P- 
sampled IC’s, on the other hand, systematically under¬ 
estimate the mass variance in spheres and the shape of 
the correlation function. Most analytic nonlinear theo¬ 
ries of structure forma tion fe.g.. lPress fc Schechteilll974 
iHamilton et all Il991i) are best understood in terms of 
real-space gravitational collapse, so it is important to 
get the real-space statistical properties correct in N- 
body simulations. Another way of making this point is 
to imagine extracting a finite size box from an underly¬ 
ing infinite density field (not worrying about periodicity, 
which comes naturally from sampling the power spec¬ 
trum). The mean density in randomly placed boxes will 
not exactly equal the cosmological mean density, and this 
sample variance should be included in TV-body simula¬ 
tions. As implied by im the most important contribu¬ 
tion of ,f-sampled IC’s is to incorporate the DC mode, so 
one could use P-sampled IC’s but treat the D C mode spe¬ 
cially using the techniques of ^ as done bv iFrenk et alJ 
lll98j^ . However, to ensure that the real-space statisti¬ 
cal properties are accurately encoded within the box, one 
should use the convolved power spectrum as described in 

SI 

Thus there are two distinct but related effects of using 
,f-sampled IC’s. To reiterate, the first is that the mean 
real-space statistical properties of an ensemble become 
unbiased. For example, for L = 50 fc,“^Mpc, ,^-sampled 
IC’s correct the 6% underestimate in as of P-sampled 
IC’s. ^-sampled IC’s also fix the shape of the correla¬ 
tion function to match the input correlation function for 
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Fig. 10.— Power spectra for the L = 100 h“^Mpc, N = 64® ensembles with P-sampled IC’s (left panels) and ^-sampled IC’s (right 
panels) for five output aun;, labeled. The black line, purple region, light green region, and red line have the same meanings as in figure!^ 
On the right side of each panel is a green line marking the expected Poisson noise contribution to the power spectrum, V/N. 


r < L/2. The second effect is that the variance of DC 
modes becomes nonzero. For example, the characteristic 
overdensity of 50 /i“^Mpc boxes is 30%! While the tests 
between the two methods in this paper have been limited 
to certain statistics, one could imagine that this latter ef¬ 
fect would be extremely important in constructing halo 
catalogs from an ensemble of simulations, for example. 

One could, of course, reduce the problem of inaccu¬ 


rately modeled real-space statistics with P-sampled IC’s 
by using very large boxes. But since each doubling of the 
box size increases computation time by about an order 
of magnitude, one could ask if it would be better to run 
one L = 100 /i“^Mpc simulation or ten realizations of 
a P = 50 /i“^Mpc ensemble. The answer is not obvi¬ 
ous and probably depends on the problem to be solved. 
The former simulation necessarily ignores the DC mode 
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variance, so erg will be biased. The latter ensemble will 
have an unbiased as but will be unable to reproduce the 
correlation function between 50 and 100 /i“^Mpc. 

I have also argued that the commonly assumed Poisson 
noise model for the discreteness of particles is not a triv¬ 
ially added contribution to the power spectrum. More 
investigation about the relationship, in the power spec¬ 
trum, between the discreteness effects and the underlying 
continuous perturbations is warranted. 

Code to generate initial conditions to sec¬ 


ond order in La grangi an perturbation the¬ 
ory l|Scoccimarrol 1 199811 is available at 
http://WWW.astro.princeton.edu/~esirko/ic. 
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ported by NASA HQ Award ^NNG04GK55G. The com¬ 
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NSF grant AST-0216105. 


APPENDIX 

DERIVATION OF THE “UNI” ^ “BOX” MAPPING 

This section outlines the derivation of the “uni” ^ “box” mapping of cosmological parameters in a A 7 ^ 0 realization 
(equations El and EJ • 

At early times, Obox = Quni and Om —> 1 , so simplifying equation 0] results in 


^box 

^box 


= H, 




O-boxI 


^box 


= H, 


0,uni 



(Al) 


and it is easily seen that (Ho.box/^^o.uni)^ = f^m,uni/f^m,box in order to accurately model the early epochs. Therefore 
let Flo.box = ^^o,uni /(1 + <(>), f^m.box = f^m,uni(l + , and OA,box = f^A.uni(l + lY where (j) and 7 are parametrizations 

of the modifications of the cosmological parameters. 

Taking the time derivative of equation El 


^box 


^^box 


^auni A/3 
ttuni 1 - A/3 


(A2) 


Inserting equations ElEl and El squaring both sides, substituting the (j) and 7 parametrizations, and then expanding 
to second order in (p, 7 , and A results in: 


n 

a 


1 -I- A -I- j + f^A(l + 27 - 2(/ -I- 7 ^ - Ajp + 
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( 1 — ‘^4’ T T ) — 
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Ha 


2 Hk 


^ ^ni I - - 
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36a2 
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a 


30^ 

2 Hk 

Hin 


5 3 

D a 

Hn 


a + 1 a^ 


Ha 




(A3) 


where 1 have dropped the “uni” subscripts on the cosmological parameters and Ouni for notational simplicity. Using 
the Eulerian form of equation m only affects some of the second-order terms. The zeroth-order equation is trivially 
satisfied. Extracting the first order equation. 


%A + Ha(27-2^)-4< 


2Hk . Hm 

3a^ 3a^ 


5 

D 



(A4) 


In order for the perturbation parameters not to have a dependence on a, 7 must equal (p. The equation then simplifies 
to 


_ 5 Hm 
“ 6 D( 1 ) 


Aq. 


(A5) 


Figiire fXTI shows the result of this approximation for Aq = 
reproduces the expected behavior for abox( 0 - 


0.4 (quite a large value!). The “box” cosmology accurately 
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Fig. A1.— Scale factor vs. cosmic time, normalized to Einstein-de Sitter expansion. The thick line shows auni vs. t for Qm = 0.27, Qa = 
0.73, h = 0.71. For this cosmology the Hubble time is I/Hq = 13.77 Gyr, and the age of the universe is 13.67 Gyr; the latter serves 
as the boundary of the plot on the right side. For Aq = 0.4, the upper dotted line is the Eulerian estimate on what a.box should be: 
abox = (1 + A)~^/^auni- The lower dotted line is the Lagrangian estimate abox = (1 ~ ^/3)auni- The thin line is Ubox vs. t for the 
perturbed cosmological parameters. 
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